Relativistic Breit–Wigner (rel_breitwigner) distribution#
The relativistic Breit–Wigner distribution is a continuous distribution on \([0,\infty)\) used in high‑energy physics to model resonances (unstable particles) via the uncertainty in their reconstructed invariant mass.
Compared to the (non‑relativistic) Breit–Wigner / Cauchy form, the relativistic version bakes in the dependence on \(x^2\) that arises naturally when modeling resonance behavior in relativistic kinematics.
Learning goals#
Understand what
rel_breitwignermodels and how it relates tocauchy/ Breit–Wigner.Work with the PDF and CDF (including a numerically robust NumPy implementation).
Compute mean/variance (and understand why higher moments diverge).
Sample from scratch (NumPy-only) via numerical inverse CDF.
Use
scipy.stats.rel_breitwignerfor evaluation, simulation, and fitting.
import platform
import numpy as np
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots
import scipy
from scipy import integrate, optimize
from scipy.stats import cauchy, rel_breitwigner
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
np.set_printoptions(precision=6, suppress=True)
rng = np.random.default_rng(7)
print("Python", platform.python_version())
print("NumPy", np.__version__)
print("SciPy", scipy.__version__)
Python 3.12.9
NumPy 1.26.2
SciPy 1.15.0
1) Title & Classification#
Name:
rel_breitwignerType: continuous distribution
Support (standardized): \(x \in [0, \infty)\)
Parameter space:
shape: \(\rho > 0\)
location: \(\text{loc} \in \mathbb{R}\)
scale: \(\text{scale} > 0\)
We write (standardized form):
SciPy uses the name rel_breitwigner and implements the usual location/scale transform:
2) Intuition & Motivation#
2.1 What it models#
In collider experiments, unstable particles (resonances) are reconstructed from their decay products. The reconstructed invariant mass is not a single number: it fluctuates due to the particle’s finite lifetime (via the uncertainty principle) and detector effects.
The relativistic Breit–Wigner distribution is a standard idealized model for this resonance line‑shape.
2.2 Typical real‑world use cases#
Resonance modeling: \(Z^0\), \(W\), \(\phi\), \(\rho\) mesons, …
Peak finding and mass/width estimation from invariant mass histograms
Likelihood-based inference in particle physics analyses
2.3 Relations to other distributions#
Cauchy / Breit–Wigner: the classic (non‑relativistic) Breit–Wigner is the Cauchy distribution.
For large \(\rho\),
rel_breitwignerlooks locally like a Cauchy near its mode (see Section 5).Heavy tails imply that not all moments exist: mean and variance exist, but skewness and kurtosis diverge.
3) Formal Definition#
SciPy’s standardized PDF is
where
CDF#
The CDF is the integral
and has a closed form that SciPy evaluates using complex arithmetic (it’s real‑valued after taking an imaginary part).
Location/scale#
With location and scale:
def _check_positive(name: str, value: float) -> float:
value = float(value)
if not np.isfinite(value) or value <= 0:
raise ValueError(f"{name} must be a finite positive number")
return value
def rel_breitwigner_k(rho: float) -> float:
"""Normalization constant k(rho) for the standardized PDF."""
rho = _check_positive("rho", rho)
return (
2
* np.sqrt(2)
* rho**2
* np.sqrt(rho**2 + 1)
/ (np.pi * np.sqrt(rho**2 + rho * np.sqrt(rho**2 + 1)))
)
def rel_breitwigner_pdf(
x: np.ndarray | float,
rho: float,
*,
loc: float = 0.0,
scale: float = 1.0,
) -> np.ndarray:
"""PDF of rel_breitwigner with optional loc/scale (NumPy-only)."""
rho = _check_positive("rho", rho)
scale = _check_positive("scale", scale)
x = np.asarray(x, dtype=float)
y = (x - float(loc)) / scale
# A numerically stable form used in SciPy: pdf(x) = (k/rho**2) / ( ((x^2-rho^2)/rho)**2 + 1 )
C = (
np.sqrt(2 * (1 + 1 / rho**2) / (1 + np.sqrt(1 + 1 / rho**2)))
* 2
/ np.pi
)
denom = (((y - rho) * (y + rho) / rho) ** 2 + 1)
pdf = (C / denom) / scale
return np.where(y >= 0, pdf, 0.0)
def rel_breitwigner_logpdf(
x: np.ndarray | float,
rho: float,
*,
loc: float = 0.0,
scale: float = 1.0,
) -> np.ndarray:
"""Log-PDF of rel_breitwigner with optional loc/scale (NumPy-only)."""
rho = _check_positive("rho", rho)
scale = _check_positive("scale", scale)
x = np.asarray(x, dtype=float)
y = (x - float(loc)) / scale
C = (
np.sqrt(2 * (1 + 1 / rho**2) / (1 + np.sqrt(1 + 1 / rho**2)))
* 2
/ np.pi
)
denom = (((y - rho) * (y + rho) / rho) ** 2 + 1)
logpdf = np.log(C) - np.log(denom) - np.log(scale)
return np.where(y >= 0, logpdf, -np.inf)
def rel_breitwigner_cdf(
x: np.ndarray | float,
rho: float,
*,
loc: float = 0.0,
scale: float = 1.0,
) -> np.ndarray:
"""CDF of rel_breitwigner with optional loc/scale (NumPy-only).
This mirrors SciPy's closed-form implementation, which uses complex arithmetic.
"""
rho = _check_positive("rho", rho)
scale = _check_positive("scale", scale)
x = np.asarray(x, dtype=float)
y = (x - float(loc)) / scale
C = np.sqrt(2 / (1 + np.sqrt(1 + 1 / rho**2))) / np.pi
y_pos = np.where(y >= 0, y, 0.0)
z = np.sqrt(-1 + 1j / rho) * np.arctan(y_pos / np.sqrt(-rho * (rho + 1j)))
cdf = C * 2 * np.imag(z)
cdf = np.clip(cdf, 0.0, 1.0)
return np.where(y >= 0, cdf, 0.0)
# Quick numerical sanity checks
rho_test = 2.0
xs = np.linspace(0, 50, 200_001)
area = np.trapz(rel_breitwigner_pdf(xs, rho_test), xs)
print("Approx integral of PDF over [0, 50]:", area)
print("CDF(0) =", rel_breitwigner_cdf(0.0, rho_test))
print("CDF(50) =", rel_breitwigner_cdf(50.0, rho_test))
Approx integral of PDF over [0, 50]: 0.9999926082589535
CDF(0) = 0.0
CDF(50) = 0.9999926082589539
4) Moments & Properties#
4.1 Mean, variance, skewness, kurtosis#
Because the tail behaves like \(f(x)\sim k/x^4\), moments exist only up to order \(n<3\):
For \(X\sim \mathrm{RelBW}(\rho)\) (standardized):
Mean \(\mathbb{E}[X]\) exists and is finite.
Variance \(\mathrm{Var}(X)\) exists and is finite.
Skewness and kurtosis are infinite (SciPy reports
nan).
4.2 MGF and characteristic function#
The MGF \(M_X(t)=\mathbb{E}[e^{tX}]\) is finite for \(t<0\) and diverges for \(t>0\) (polynomial tails cannot compete with \(e^{tX}\)).
The characteristic function \(\varphi_X(t)=\mathbb{E}[e^{itX}]\) always exists.
4.3 Entropy#
The differential entropy
does not have a commonly used simple closed form; SciPy evaluates it numerically.
def rel_breitwigner_mean(rho: float) -> float:
rho = _check_positive("rho", rho)
A = np.sqrt(rho**2 + 1)
return (A / np.pi) * np.sqrt(2 * rho / (rho + A)) * (np.pi / 2 + np.arctan(rho))
def rel_breitwigner_second_moment(rho: float) -> float:
rho = _check_positive("rho", rho)
return rho * np.sqrt(rho**2 + 1)
def rel_breitwigner_var(rho: float) -> float:
m = rel_breitwigner_mean(rho)
return rel_breitwigner_second_moment(rho) - m**2
for rho in [0.2, 1.0, 3.5]:
m_formula = rel_breitwigner_mean(rho)
v_formula = rel_breitwigner_var(rho)
m_scipy, v_scipy = rel_breitwigner.stats(rho, moments="mv")
print(f"rho={rho}")
print(" mean (formula, SciPy):", m_formula, m_scipy)
print(" var (formula, SciPy):", v_formula, v_scipy)
print("entropy(rho=1.0) via SciPy:", rel_breitwigner.entropy(1.0))
rho=0.2
mean (formula, SciPy): 0.3286859777676976 0.3286859777676975
var (formula, SciPy): 0.09592630856260402 0.09592630856260406
rho=1.0
mean (formula, SciPy): 0.9653913793583742 0.965391379358374
var (formula, SciPy): 0.4822330470336309 0.4822330470336309
rho=3.5
mean (formula, SciPy): 3.284899599858133 3.2848995998581327
var (formula, SciPy): 1.9496269250927831 1.9496269250927867
entropy(rho=1.0) via SciPy: 0.7677961497524766
def rel_breitwigner_cf(t: float, rho: float) -> complex:
"""Characteristic function via numerical integration (SciPy quad)."""
rho = _check_positive("rho", rho)
t = float(t)
f_cos = lambda x: np.cos(t * x) * rel_breitwigner_pdf(x, rho)
f_sin = lambda x: np.sin(t * x) * rel_breitwigner_pdf(x, rho)
re, _ = integrate.quad(f_cos, 0, np.inf, limit=300)
im, _ = integrate.quad(f_sin, 0, np.inf, limit=300)
return re + 1j * im
for t in [0.0, 0.5, 1.0, 2.0]:
val = rel_breitwigner_cf(t, rho=2.0)
print(f"phi({t}) = {val}")
phi(0.0) = (0.9999999999999997+0j)
phi(0.5) = (0.5630438239678969+0.7174512887274822j)
phi(1.0) = (-0.1597740266239293+0.7104281135014987j)
phi(2.0) = (-0.2863681210052288-0.20255826100750324j)
5) Parameter Interpretation#
5.1 Meaning of parameters#
In the physics motivation, for a resonance with characteristic mass \(M_0\) and decay width \(\Gamma\) (in natural units), SciPy’s parametrization uses
shape: \(\rho = M_0/\Gamma\)
scale: \(\text{scale}=\Gamma\)
location: typically \(\text{loc}=0\)
so that the mode occurs at
5.2 How the shape changes#
Increasing \(\rho\) moves the peak to the right and makes the distribution more “Cauchy‑like” near its mode.
Increasing scale stretches the distribution (wider peak; larger absolute variance).
def plot_pdf_cdf_for_rhos(rhos, x_max=10.0):
xs = np.linspace(0, x_max, 900)
fig = make_subplots(
rows=1,
cols=2,
subplot_titles=("PDF", "CDF"),
horizontal_spacing=0.12,
)
for rho in rhos:
fig.add_trace(
go.Scatter(
x=xs,
y=rel_breitwigner_pdf(xs, rho),
mode="lines",
name=f"rho={rho}",
),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(
x=xs,
y=rel_breitwigner_cdf(xs, rho),
mode="lines",
name=f"rho={rho}",
showlegend=False,
),
row=1,
col=2,
)
fig.update_xaxes(title_text="x", row=1, col=1)
fig.update_xaxes(title_text="x", row=1, col=2)
fig.update_yaxes(title_text="f(x)", row=1, col=1)
fig.update_yaxes(title_text="F(x)", row=1, col=2)
fig.update_layout(title="rel_breitwigner: PDF and CDF (standardized)")
return fig
plot_pdf_cdf_for_rhos([0.5, 1.0, 2.0, 5.0], x_max=12)
# Local Cauchy approximation for large rho near the mode.
# For large rho, around x pprox rho, rel_breitwigner(x; rho) pprox Cauchy(loc=rho, scale=1/2).
rho = 10.0
xs = np.linspace(max(0, rho - 3), rho + 3, 800)
pdf_rbw = rel_breitwigner_pdf(xs, rho)
pdf_cauchy = cauchy.pdf(xs, loc=rho, scale=0.5)
fig = go.Figure()
fig.add_trace(go.Scatter(x=xs, y=pdf_rbw, mode="lines", name="rel_breitwigner"))
fig.add_trace(
go.Scatter(
x=xs,
y=pdf_cauchy,
mode="lines",
name="Cauchy approx",
line=dict(dash="dash"),
)
)
fig.update_layout(
title="Large-rho local approximation near the mode",
xaxis_title="x",
yaxis_title="density",
)
fig
6) Derivations#
6.1 Expectation (mean)#
Start from
Substitute \(u=x^2\) so that \(du=2x\,dx\):
This is a shifted Cauchy kernel. Using
with \(a=\rho^2\) and \(b=\rho\) gives
6.2 Second moment and variance#
The second raw moment is
This integral can be evaluated in closed form (e.g. by extending to an even integral over \(\mathbb{R}\) and applying residue calculus), yielding a remarkably simple result:
Therefore
Higher moments diverge because \(f(x)\sim k/x^4\).
6.3 Likelihood#
Given i.i.d. samples \(x_1,\dots,x_n\) from the standardized model (support \(x_i\ge 0\)), the log-likelihood for \(\rho\) is
With location/scale, let \(y_i=(x_i-\text{loc})/\text{scale}\) and add the Jacobian term \(-n\log(\text{scale})\).
def nll_rho_only(rho: float, data: np.ndarray) -> float:
"""Negative log-likelihood for rho (standardized: loc=0, scale=1)."""
rho = float(rho)
if rho <= 0:
return np.inf
x = np.asarray(data, dtype=float)
if np.any(x < 0):
return np.inf
k = rel_breitwigner_k(rho)
ll = len(x) * np.log(k) - np.sum(np.log((x**2 - rho**2) ** 2 + rho**2))
return -ll
# Synthetic MLE demo
rho_true = 2.5
x = rel_breitwigner.rvs(rho_true, size=5_000, random_state=rng)
res = optimize.minimize_scalar(
nll_rho_only,
bounds=(1e-3, 50.0),
args=(x,),
method="bounded",
)
print("true rho:", rho_true)
print("MLE rho :", res.x)
true rho: 2.5
MLE rho : 2.4911323685576505
7) Sampling & Simulation (NumPy-only)#
SciPy can sample from rel_breitwigner, but here we implement a from-scratch sampler that uses only:
a source of uniforms \(U\sim\mathrm{Uniform}(0,1)\) (NumPy RNG)
the CDF \(F(x\mid\rho)\)
numerical inversion (bisection)
Algorithm (inverse CDF via bisection)#
For each \(u\in(0,1)\):
Find a bracket \([\ell, h]\) with \(F(\ell)\le u \le F(h)\) (start at \(\ell=0\) and expand \(h\) by doubling).
Repeat bisection for a fixed number of iterations.
This is robust and easy to implement for any monotone CDF.
def rel_breitwigner_ppf_bisect(
u: np.ndarray | float,
rho: float,
*,
tol: float = 1e-12,
max_iter: int = 80,
) -> np.ndarray:
"""Inverse CDF (PPF) via bisection (NumPy-only)."""
rho = _check_positive("rho", rho)
u = np.asarray(u, dtype=float)
if np.any(~np.isfinite(u)):
raise ValueError("u must be finite")
out = np.empty_like(u)
mask0 = u <= 0
mask1 = u >= 1
mask = ~(mask0 | mask1)
out[mask0] = 0.0
out[mask1] = np.inf
if not np.any(mask):
return out
uu = u[mask]
lo = np.zeros_like(uu)
hi = np.full_like(uu, max(1.0, rho))
# Expand hi until all targets are bracketed
for _ in range(100):
cdf_hi = rel_breitwigner_cdf(hi, rho)
need = cdf_hi < uu
if not np.any(need):
break
hi[need] *= 2
else:
raise RuntimeError("Failed to bracket all u values")
# Bisection
for _ in range(max_iter):
mid = 0.5 * (lo + hi)
cdf_mid = rel_breitwigner_cdf(mid, rho)
go_right = cdf_mid < uu
lo[go_right] = mid[go_right]
hi[~go_right] = mid[~go_right]
if np.max(hi - lo) < tol:
break
out[mask] = 0.5 * (lo + hi)
return out
def rel_breitwigner_rvs_numpy(rho: float, size: int, *, rng: np.random.Generator) -> np.ndarray:
u = rng.random(size)
return rel_breitwigner_ppf_bisect(u, rho)
# Quick check: sample moments vs theory
rho = 2.0
x_samp = rel_breitwigner_rvs_numpy(rho, 30_000, rng=rng)
print("sample mean :", x_samp.mean())
print("theory mean :", rel_breitwigner_mean(rho))
print("sample var :", x_samp.var())
print("theory var :", rel_breitwigner_var(rho))
sample mean : 1.8575346488861357
theory mean : 1.8521891046099694
sample var : 0.8987076883363847
theory var : 1.041531475763699
8) Visualization#
We’ll visualize:
the PDF and CDF for different \(\rho\)
a Monte Carlo histogram compared to the theoretical PDF
rho = 2.0
xs = np.linspace(0, 15, 1200)
samples = rel_breitwigner_rvs_numpy(rho, 50_000, rng=rng)
fig = make_subplots(
rows=1,
cols=2,
subplot_titles=("PDF + Monte Carlo histogram", "CDF"),
horizontal_spacing=0.12,
)
fig.add_trace(
go.Histogram(
x=samples,
nbinsx=120,
histnorm="probability density",
name="samples",
opacity=0.55,
),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(
x=xs,
y=rel_breitwigner_pdf(xs, rho),
mode="lines",
name="theoretical pdf",
line=dict(width=2),
),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(
x=xs,
y=rel_breitwigner_cdf(xs, rho),
mode="lines",
name="cdf",
),
row=1,
col=2,
)
fig.update_xaxes(title_text="x", row=1, col=1)
fig.update_xaxes(title_text="x", row=1, col=2)
fig.update_yaxes(title_text="density", row=1, col=1)
fig.update_yaxes(title_text="probability", row=1, col=2)
fig.update_layout(title=f"rel_breitwigner (rho={rho})")
fig
9) SciPy Integration (scipy.stats.rel_breitwigner)#
SciPy exposes the distribution as scipy.stats.rel_breitwigner.
Key methods:
pdf,logpdfcdf,ppf(numerical inversion)rvsfit
Physical tip: when fitting resonance masses, it is usually appropriate to fix loc=0 (see SciPy docs).
rho = 2.0
x0 = np.array([0.0, 0.5, 1.0, 2.0, 5.0])
print("pdf:", rel_breitwigner.pdf(x0, rho))
print("cdf:", rel_breitwigner.cdf(x0, rho))
s = rel_breitwigner.rvs(rho, size=5, random_state=rng)
print("rvs:", s)
# Fit demo (standardized, so loc should be ~0 and scale ~1)
data = rel_breitwigner.rvs(2.5, size=5_000, random_state=rng)
rho_hat, loc_hat, scale_hat = rel_breitwigner.fit(data, floc=0)
print("fit rho, loc, scale:", rho_hat, loc_hat, scale_hat)
pdf: [0.138329 0.153167 0.212814 0.691646 0.006217]
cdf: [0. 0.071569 0.160358 0.583635 0.99095 ]
rvs: [2.292835 3.588495 0.863896 2.080456 0.086615]
fit rho, loc, scale: 2.5159002307382545 0 1.0010439852719006
# Example: Z0 resonance parameters (Particle Data Group numbers in SciPy docs)
M0 = 91.1876
Gamma = 2.4952
rho = M0 / Gamma
xs = np.linspace(70, 110, 1200)
pdf = rel_breitwigner.pdf(xs, rho, loc=0, scale=Gamma)
fig = go.Figure()
fig.add_trace(go.Scatter(x=xs, y=pdf, mode="lines", name="Z0 model"))
fig.update_layout(
title="Relativistic Breit–Wigner example (Z0 line shape)",
xaxis_title="invariant mass M (GeV)",
yaxis_title="density",
)
fig
10) Statistical Use Cases#
10.1 Hypothesis testing (signal vs background)#
A common pattern is testing for the presence of a resonant signal on top of a smoother background.
We’ll do a toy likelihood‑ratio test between:
\(H_0\): background only (uniform)
\(H_1\): mixture of background +
rel_breitwignersignal (known parameters)
10.2 Bayesian modeling#
You can treat rel_breitwigner as a likelihood for \(\rho\) (and possibly scale), assign priors, and compute a posterior.
10.3 Generative modeling#
In simulation pipelines, rel_breitwigner acts as a parametric generator of resonance masses.
def logpdf_uniform_bg(x: np.ndarray, L: float) -> np.ndarray:
x = np.asarray(x, dtype=float)
return np.where((0 <= x) & (x <= L), -np.log(L), -np.inf)
def logpdf_signal(x: np.ndarray, rho: float, scale: float) -> np.ndarray:
return rel_breitwigner_logpdf(x, rho, loc=0.0, scale=scale)
def logpdf_mixture(x: np.ndarray, w: float, rho: float, scale: float, L: float) -> np.ndarray:
# log( w*s + (1-w)*b ) computed stably
x = np.asarray(x, dtype=float)
log_s = logpdf_signal(x, rho, scale)
log_b = logpdf_uniform_bg(x, L)
a = np.log(w) + log_s
b = np.log(1 - w) + log_b
m = np.maximum(a, b)
return m + np.log(np.exp(a - m) + np.exp(b - m))
# Toy LRT
rng = np.random.default_rng(7)
L = 200.0
rho = 36.5
scale = 2.5
w = 0.15
n = 400
# observed data under H1
x_sig = rel_breitwigner.rvs(rho, loc=0, scale=scale, size=int(n * w), random_state=rng)
x_bg = rng.uniform(0, L, size=n - len(x_sig))
x_obs = np.concatenate([x_sig, x_bg])
ll0 = np.sum(logpdf_uniform_bg(x_obs, L))
ll1 = np.sum(logpdf_mixture(x_obs, w=w, rho=rho, scale=scale, L=L))
T_obs = 2 * (ll1 - ll0)
print("Observed LRT statistic T:", T_obs)
# Monte Carlo null distribution under H0
B = 400
T_null = np.empty(B)
for b in range(B):
x0 = rng.uniform(0, L, size=n)
ll0_b = np.sum(logpdf_uniform_bg(x0, L))
ll1_b = np.sum(logpdf_mixture(x0, w=w, rho=rho, scale=scale, L=L))
T_null[b] = 2 * (ll1_b - ll0_b)
p_value = np.mean(T_null >= T_obs)
print("Monte Carlo p-value:", p_value)
fig = go.Figure()
fig.add_trace(go.Histogram(x=T_null, nbinsx=40, name="T under H0"))
fig.add_vline(x=T_obs, line_width=2, line_dash="dash", line_color="black")
fig.update_layout(
title="Toy likelihood-ratio test: background vs (background+signal)",
xaxis_title="T",
yaxis_title="count",
)
fig
Observed LRT statistic T: 96.15435313231501
Monte Carlo p-value: 0.0
# Simple Bayesian inference for rho (standardized: scale=1, loc=0)
# Posterior p(rho | x) ∝ p(x | rho) p(rho)
rng = np.random.default_rng(7)
rho_true = 2.0
x = rel_breitwigner.rvs(rho_true, size=400, random_state=rng)
rho_grid = np.linspace(0.2, 8.0, 600)
# Prior: log-normal on rho (weakly informative)
mu, sigma = 0.0, 0.7
log_prior = -np.log(rho_grid * sigma * np.sqrt(2 * np.pi)) - (np.log(rho_grid) - mu) ** 2 / (2 * sigma**2)
log_like = np.array([-nll_rho_only(r, x) for r in rho_grid])
log_post_unnorm = log_like + log_prior
log_post = log_post_unnorm - np.max(log_post_unnorm)
post = np.exp(log_post)
post = post / np.trapz(post, rho_grid)
rho_map = rho_grid[np.argmax(post)]
print("true rho:", rho_true)
print("MAP rho :", rho_map)
fig = go.Figure()
fig.add_trace(go.Scatter(x=rho_grid, y=post, mode="lines", name="posterior"))
fig.add_vline(x=rho_true, line_width=2, line_dash="dash", line_color="black")
fig.add_vline(x=rho_map, line_width=2, line_dash="dot", line_color="gray")
fig.update_layout(
title="Posterior over rho (scale fixed to 1)",
xaxis_title="rho",
yaxis_title="density",
)
fig
true rho: 2.0
MAP rho : 2.0230383973288815
# Generative modeling: draw synthetic resonance masses (Z0 example)
rng = np.random.default_rng(7)
M0 = 91.1876
Gamma = 2.4952
rho = M0 / Gamma
masses = Gamma * rel_breitwigner_rvs_numpy(rho, 30_000, rng=rng)
fig = go.Figure()
fig.add_trace(
go.Histogram(
x=masses,
nbinsx=140,
name="simulated events",
histnorm="probability density",
opacity=0.7,
)
)
xs = np.linspace(70, 110, 900)
fig.add_trace(
go.Scatter(
x=xs,
y=rel_breitwigner.pdf(xs, rho, loc=0, scale=Gamma),
mode="lines",
name="model pdf",
)
)
fig.update_layout(
title="Generative simulation: Z0 resonance masses",
xaxis_title="invariant mass M (GeV)",
yaxis_title="density",
)
fig
11) Pitfalls#
Invalid parameters:
rhoandscalemust be strictly positive.Support matters: in standardized form, \(x<0\) has density 0; with location, support is \(x\ge\text{loc}\).
Higher moments diverge: sample skewness/kurtosis can be unstable or meaningless.
Fitting
loc: allowinglocto vary can produce non-physical negative mass support; in resonance modeling, preferfloc=0.Numerical CDF/PPF: SciPy’s
ppfrelies on numerical inversion; for extreme probabilities it can be slow. For simulation,rvsis usually preferable.
12) Summary#
rel_breitwigneris a continuous distribution on \([0,\infty)\) used to model relativistic resonance line shapes.PDF: \(f(x\mid\rho)=\dfrac{k(\rho)}{(x^2-\rho^2)^2+\rho^2}\) with \(\rho>0\).
Mean and variance exist, but skewness/kurtosis diverge.
Sampling can be done via numerical inverse CDF (NumPy-only) or via SciPy’s
rvs.In physics parametrization: \(\rho=M_0/\Gamma\) and
scale=Γ, so the mode is at \(M_0\).